Get blusky raster for all days

bluesky_grid <- bluesky_load(file = blueskyPath)

bluesky_oct09a <- createBlueskyRasterSlice(bluesky_grid, 2017100921)
bluesky_oct09t <- createBlueskyRasterSlice(bluesky_grid, 2017100919)
bluesky_oct10a <- createBlueskyRasterSlice(bluesky_grid, 2017101021)
bluesky_oct10t <- createBlueskyRasterSlice(bluesky_grid, 2017101019)
bluesky_oct11a <- createBlueskyRasterSlice(bluesky_grid, 2017101121)
bluesky_oct11t <- createBlueskyRasterSlice(bluesky_grid, 2017101119)
bluesky_oct12a <- createBlueskyRasterSlice(bluesky_grid, 2017101221)
bluesky_oct12t <- createBlueskyRasterSlice(bluesky_grid, 2017101219)
bluesky_oct13a <- createBlueskyRasterSlice(bluesky_grid, 2017101321)
bluesky_oct13t <- createBlueskyRasterSlice(bluesky_grid, 2017101319)

Resample maiac data to match bluesky

rasterList <- list(AQUA_09, TERRA_09, AQUA_10, TERRA_10, AQUA_11, TERRA_11, AQUA_12, TERRA_12, AQUA_13, TERRA_13)
rslist <- list()
for ( i in 1:5 ) {
  AQUA <- rasterList[[i*2 - 1]]$corrected
  TERRA <- rasterList[[i*2]]$corrected
  
  # Get the day of the month as a zero-padded string for printing and selecting the correct value from monitoring data
  day <- stringr::str_pad(i+8, 2, "left", "0")
  
  # Resample aqua and terra
  aqua_r <- resample(AQUA, bluesky_oct09a)
  assign(paste0("rs_", day, "a"), aqua_r)
  
  terra_r <- resample(TERRA, bluesky_oct09t)
  assign(paste0("rs_", day, "t"), terra_r)
  
  rslist[[paste0("rs_", day, "a")]] <- aqua_r
  rslist[[paste0("rs_", day, "t")]] <- terra_r
}

Compare bluesky (converted to AOD) and MAIAC AOD values

blueskyList <- list(bluesky_oct09a, bluesky_oct09t, bluesky_oct10a, bluesky_oct10t, bluesky_oct11a, bluesky_oct11t,
                    bluesky_oct12a, bluesky_oct12t, bluesky_oct13a, bluesky_oct13t)

par(pty = "s")
# Separate plots for each day
for ( i in 1:5 ) {
  AQUA <- rslist[[i*2 - 1]]$corrected
  TERRA <- rslist[[i*2]]$corrected
  bsA <- blueskyList[[i*2 - 1]]
  bsT <- blueskyList[[i*2]]
  
  # Get the day of the month as a zero-padded string for printing and selecting the correct value from monitoring data
  day <- stringr::str_pad(i+8, 2, "left", "0")
  
  # plot Aqua values
  plot(toAOD(values(bsA))~values(AQUA), xlab = "MAIAC", ylab = "Bluesky", main = paste0("Aqua Oct ", day), xlim = c(0, 3), ylim = c(0, 3))

  
  # plot Terra values
  plot(toAOD(values(bsT))~values(TERRA), xlab = "MAIAC", ylab = "Bluesky", main = paste0("Terra Oct ", day), xlim = c(0, 3), ylim = c(0, 3))
}

Make some side-by-side maps

# automatically-generated color palatte
colors <- colorRampPalette(c("white", "orange", "red"))(9)

# define breaks for colors
breaks <- c(0, 0.03, 0.1, 0.2, 0.3, 0.5, 0.8, 4)

difbreaks <- c(-5, -2, -1, -.8, -.5, -.2, -.1, .1, .2, .5, .8, 1, 2, 10)
difcolors <- colorRampPalette(c("blue", "white", "red"))(13)

par(mfcol = c(1,3))

for ( i in 1:5 ) {
  AQUA <- rslist[[i*2 - 1]]$corrected
  TERRA <- rslist[[i*2]]$corrected
  bsA <- blueskyList[[i*2 - 1]]
  values(bsA) <- toAOD(values(bsA))
  bsT <- blueskyList[[i*2]]
  values(bsT) <- toAOD(values(bsT))
  
  # Get the day of the month as a zero-padded string for printing
  day <- stringr::str_pad(i+8, 2, "left", "0")
  
  # Map bluesky image
  plot(bsT, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("bluesky Oct ", day, ", 2017"), legend = FALSE)
  plot(USCensusStates, add = TRUE)
  
  # Map satellite image
  plot(TERRA, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("TERRA ", "Oct ", day, ", 2017"), legend = FALSE)
  plot(USCensusStates, add = TRUE)
  
  # Map the difference
  plot(bsT-TERRA, breaks = difbreaks, col = difcolors, colNA = "gray90", main = "Difference in AOD bluesky to MAIAC", legend = FALSE)
  plot(USCensusStates, add = TRUE)
  
  
  # Do the same for AQUA
   plot(bsA, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("bluesky Oct ", day, ", 2017"), legend = FALSE)
  plot(USCensusStates, add = TRUE)
  plot(AQUA, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("AQUA ", "Oct ", day, ", 2017"), legend = FALSE)
  plot(USCensusStates, add = TRUE)
  plot(bsA-AQUA, breaks = difbreaks, col = difcolors, colNA = "gray90", main = "Difference AOD bluesky to MAIAC", legend = FALSE)
  plot(USCensusStates, add = TRUE)
  
}